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The computation of the entire Lyapunov spectrum for ex- 
tended dynamical systems is a very time consuming task. If 
the system is in a chaotic spatio-temporal regime it is pos- 
sible to approximately reconstruct the Lyapunov spectrum 
from the spectrum of a sub-system in a very cost effective 
way. In this work we present a new rescaling method, which 
gives a significantly better fit to the original Lyapunov spec- 
trum. It is inspired by the stability analysis of the homo- 
geneous evolution in a one-dimensional coupled map lattice 
but appears to be equally valid in a much wider range of 
cases. We evaluate the performance of our rescaling method 
by comparing it to the conventional rescaling (dividing by 
the relative sub-system volume) for one and two-dimensional 
lattices in spatio-temporal chaotic regimes. In doing so we 
notice that the Lyapunov spectra for consecutive sub-system 
sizes are interleaved and we discuss the possible ways in 
which this may arise. Finally, we use the new rescaling 
to approximate quantities derived from the Lyapunov spec- 
trum (largest Lyapunov exponent, Lyapunov dimension and 
Kolmogorov- Sinai entropy) finding better convergence as the 
sub-system size is increased than with conventional rescal- 
ing. 



CONTENTS 



Introduction 



[I 


One-dimensional coupled map lattices 


3 


|A Interleaving and rescaling for homoge- 








neons states 




3 




P Interleaving and rescaling for couplccj 








logistic maps 




6 




|C Estimation of quantities derived from 








the Lyapunov spectrum 




10 



III More general extended dynamical sys- 



terns 



A Chaotic neural networks 



B Two-dimensional logistic lattice 



G Host-parasitoid system 



[V Conclusions 



11 

11 

12 
14 

16 



I. INTRODUCTION 



Spatio-temporal systems give rise to a wide range of 
interesting phenomena that cannot occur in dynami- 
cal systems with only a few degrees of freedom. The 
most common approach to modelling complex spatio- 
temporal behaviour is through the use of partial differ- 
ential equations (PDE's). The analysis and even the 
numerical integration of PDE's is usually quite intri- 
cate. Thus, if one desires to study the full range of com- 
plex spatio-temporal behaviour whilst conserving a rel- 
atively simple dynamical framework, a better approach 
is to consider discrete spatio-temporal systems. By this 
we mean a collection of coupled simple low-dimensional 
dynamical units arranged on a spatial lattice. The cou- 
pling is usually (but not always) restricted to a finite 
neighbourhood. An immediate advantage of such sys- 
tems is their straightforward computational implemen- 
tation. Another possible advantage is that the local dy- 
namics at each lattice site in the uncoupled limit can be 
thoroughly analysed. The knowledge of such local dy- 
namics in the uncoupled limit can help to provide some 
insight of the complexity of the coupled system. 

In this paper we are particularly interested in the char- 
acterization of chaos in such extended dynamical sys- 
tems. The most basic tool for analyzing a chaotic sys- 
tem is its Lyapunov exponents. The Lyapunov expo- 
nents are an important invariant of nonlinear dynami- 
cal systems and are closely related to other quantities 
of interest. Consider a discrete spatio-temporal system 
with N state variables where N is the number of lo- 
cal variables times the spatial volume of the system, for 
example N — r\L d for a d-dimensional cubic lattice of 
side L with rj variables in each node. For such an N- 
dimensional system there exist N Lyapunov exponents 
corresponding to the rates of expansion and/or contrac- 
tion of nearby orbits in the tangent space in each di- 
mension. The Lyapunov spectrum (LS) is defined as the 
set °f the N Lyapunov exponents arranged in 

decreasing order. The LS is very useful in the character- 
ization of a chaotic attractor since it gives an estimate 
of its dimension by means of the Lyapunov dimension 
Dl (Kaplan- Yorke conjecture |l|,||) defined as 
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D L =j + ^Y J \ i , (1) 



where j is the largest integer for which Yji=i > 0- 
Another useful invariant that can be derived from the 
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LS is the so called Kolmogorov-Sinai (KS) entropy h 
that can be bounded from above by the sum of the pos- 
itive Lyapunov exponents and that in many cases 
can be well approximated by || 



(2) 



The KS entropy quantifies the mean rate of information 
production in a system, or alternatively the mean rate 
of growth of uncertainty in a system subjected to small 
perturbations. 

When dealing with extended dynamical systems, the 
high number of variables, and even the number of ef- 
fective degrees of freedom, often leads to severe difficul- 
ties because of the large amount of resources (comput- 
ing time and memory space) required for many com- 
putations. Therefore it is useful, and often crucial, to 
develop techniques that derive information about the 
whole system by analyzing a comparatively small sub- 
system. For dynamical systems with only a few degrees 
of freedom the computation of the LS is a straightfor- 
ward task; however, when the number of degrees of 
freedom gets large {e.g. a few hundred) it becomes a 
painstaking process In particular, any algorithm 

to compute the LS must contain two fundamental proce- 
dures; one to multiply by the Jacobian at each time step 
and the other to perform some kind of reorthonormal- 
ization j?J. The latter is required to prevent the Jaco- 
bian matrix progressively getting more ill-conditioned, 
until the largest Lyapunov exponent swamps all the 
others. Such orthogonalization procedures are based 
upon the factorization of the Jacobian matrix into a 
product of an orthogonal matrix Q and an upper tri- 
angular matrix R. The two most widespread meth- 
ods for achieving such orthogonalization are based upon 
modified Gram-Schmidt (MGS) orthogonalization and 
the so-called HQR decomposition that uses Householder 
transformations. The MGS-based methods are widely 
used because of their quite simple numerical implemen- 
tation though they are known to introduce small errors 
due to the fact that the orthogonality of the matrix Q 
may fail. The HQR-based methods are more difficult to 
implement but they give a better approximation of the 
LS H since they do not have the problem of losing or- 
thogonality of the matrix Q. The difficulty in using any 
of these methods for computing the LS of systems with 
a high number of degrees of freedom N is that they re- 
quire 0(N 3 ) operations ||. The usual naive algorithm 
for matrix multiplication is also 0(N 3 ), so that overall 
computing the full LS is an 0(N 3 ) process (in principle 
matrix multiplication can be done faster than 0(N 3 ) us- 
ing specialized techniques, but this hardly seems worth 
doing under the circumstances). As an example, the 
computation of the LS using a HQR method for a logis- 
tic coupled map lattice with N — 100 takes a few hours 
on a standard workstation. When the system size is an 
order of magnitude larger (e.g. for two or more spatial 
dimensions) and/or the convergence of the Lyapunov 
exponents is rather slow, the task quickly becomes in- 
feasible. Therefore one must rely on other techniques 



to approximate the LS for large systems. 

One such technique to estimate the LS in a fully spatio- 
temporal chaotic regime is to take a principal sub- 
matrix of the Jacobian and compute the LS for this 
sub-system. It has been observed in a wide range of 
spatio-temporal systems that such a sub-system LS con- 
verges to the spectrum of the whole system under ap- 
propriate rescaling. In a number of specialized cases 
e.g. turbulent Navier-Stokes flows || and hard sphere 
gases there are rigorous results for this phe- 

nomenon. However is seems difficult to prove its oc- 
currence more generally, and certainly there are many 
systems where it is observed numerically but no rig- 
orous analysis exists. These include coupled logistic 
maps |fj|, chaotic neural networks fl2]| , coupled map lat- 
tices Jl3|,[l4| , reaction-diffusion systems jl5],[l6) (lattice of 
ODEs), turbulent fluids (fjj], the Kuramoto-Sivashinsky 
model Q (PDE's), and others. 

Such a rescaling approach consists of evolving the whole 
TV-dimensional system under the equations of motion, 
taking a sub-system of size N s to compute the LS and 
then rescaling it to obtain an estimate of the original LS. 
This method relies on the linear increase of Lyapunov 
dimension Dl and KS entropy h with the sub-system 
size (see above references) . A physical interpretation of 
this phenomenon can be given in terms of the thermody- 
namic limit of the system. A spatio-temporal system in 
a fully chaotic regime will possess a typical correlation 
length £ such that elements further apart than £ evolve 
almost independently from each other. The whole sys- 
tem can then be thought of in some sense as the union 
of several almost independent sub-systems of size £. In 
the limit when these sub-systems are completely uncou- 
pled the LS repeats itself in each one of them. If an 
interaction between the sub-systems is introduced, one 
may expect the overall LS not to be significantly al- 
tered. Thus in the limit of a large number of degrees of 
freedom, a number of Lyapunov exponents per £- volume 
may be defined. One expects such an intuitive picture 
to become more accurate in the limit of a large number 
of degrees of freedom and a small correlation length. 

When examining closer the Lyapunov spectra in the 
fully chaotic regime for several spatio-temporal systems 
we found that the Lyapunov exponents of two consec- 
utive sub-system sizes N s and N s + 1 were interleaved. 
In other words, the zth Lyapunov exponent for the sub- 
system N s lies between the ith and (i + l)th Lyapunov 
exponents of the sub-system N s + 1. The interleaving of 
the eigenvalues for a single matrix is a well-known fact 
(Cauchy's interlace theorem) and is common in many 
areas such as Sturm sequences of polynomials |L9| . Un- 
fortunately there appears to be no obvious generaliza- 
tion which would imply the same fact for sub-system 
Lyapunov Spectra. 

This paper is organized as follows. In order to study 
interleaving we begin by examining the properties of 
sub-sys tem LS in coupled map lattices in section ||. In 
section [I A we investigate the simplest case of homoge- 
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neous evolution where we are able to prove rigorously 
that interleaving and rescaling occur. This example also 
suggest a different rescaling of the sub-system LS that is 
supe rior to that hitherto used in the literature. In sec- 
tion [I B we study the interleaving and rescaling prop- 



erties of the sub-system LS in the fully chaotic regime 
for a coupled map lattice. Applying the new rescaling 
obtained from the homogeneous case turns out to lead 
to a much better fit to the whole LS. In section II C 



we show that by using this rescaling it is possible to 
extrapolate the whole LS and extract better estimates 
for the largest Lyapunov exponent, the Lyapunov di- 
mension and the KS entropy. In section ^ we examine 
the interleaving and rescaling for more complex spatio- 
temporal systems. We notice that the interleaving is 
not always exact but the proportion of Lyapunov ex- 
ponents that do not interleave is very small. We also 
present some results for two-dimensional systems and 
point out that one needs to be careful about the choice 
of sub-system variables. Finally, in section |iy| we give 
a brief recapitulation of the results together with a dis- 
cussion on the applicability of interleaving and rescaling 
to more general extended dynamical systems. 



ONE-DIMENSIONAL COUPLED MAP 
LATTICES 



II. 



Coupled map lattices |2Jj,|2l|] (CML's) are a popular 
choice for the study of fully-developed turbulence and 
pattern formation. The appeal of CML's is due on one 
hand to their computational simplicity and on the other 
to the fact that they display a wide variety of spatio- 
temporal phenomena ranging from spatio-temporal pe- 
riodic states p2|, |23 | and travelling interfaces p4 25 to 
intermittency pdfland turbulence p^J23]. A CML is a 



discrete space-time dynamical system with a continuous 
state space, in contrast to cellular automata where the 
state space is discrete. Let us denote by x? the state 
of the ith site at time n, where the integer index i runs 
from 1 to N. The CML dynamics is defined by 



(l-s)f(x?)+ e k fW +k ), 



(3) 



where we use periodic boundary conditions, / is a real 
function and we ask ^2 £fc = e as a conservation law. 
The general CML (||) couples / > left neighbours and 
r > right neighbours with coefficients 



Interleaving and rescaling for homogeneous 
states 



In order to gain some insight into interleaving and 
rescaling behaviour of the Lyapunov spectrum in ex- 
tended dynamical systems let us start with the simplest 



case of all: homogeneous evolution. We define homoge- 
neous states as states of the form X n — {a;"}^ where 
x\ L = x n is the same for all i. It is trivial that by setting 
the initial state of the lattice to a homogeneous state 
x® — x° one has that X n = {f n (x )} for all i at any 
future time n. In other words the homogeneity of the 
initial state is preserved under iteration by ([}]). 

Let us take a simple form for the coupling by using the 
most widespread model of a CML, the so called diffusive 
CML: 

x" +1 = (1 - s)f(x-) + I (f(xU) + , (4) 

where now the coupling is symmetric and only between 
nearest neighbours. We shall perform a linear stability 
analysis of homogeneous states in this system. Such an 
analysis for more general CML's has also served as the 
starting point for the study of signal propagation ]2l| 
and pattern formation p~2[ . Since (||) preserves homo- 
geneity under iteration it is natural to ask whether the 
stability of / completely determines the stability of the 
homogeneous state. The answer turns out to be yes. 

The Lyapunov exponents are given by the logarithms 
of the eigenvalues of the matrix 

T= lim [P{n) u ' ■ P(n)] 1/2n 



(5) 



where 



P(n) = J(n) ■ J(n -!)••■ J(2) • J(l) 



and where J(s) is the Jacobian matrix of the CML dy- 
namics at time s and ( • ) tr denotes matrix transpose. 
The existence of the limit in equation (|^) for almost ev- 
ery orbit (with respect to an ergodic invariant measure) 
is guaranteed by the multiplicative ergodic theorem |3(J . 
For the homogeneous lattice 

J(n) = fj. n - M (6) 

where fi n = f'(x n ) is the multiplier of the local map 
and M is the constant matrix 



M 



fl-E 

e/2 


V e/2 



e/2 
1-e e/2 
e/2 1 - e 



e/2 \ 





e/2 1 -e J 



The matrix M is not only symmetric but also circulant. 
Recall that a matrix is circulant if in each successive 
row the elements move to the right one position (with 
wrap-around at the edges) Q . It is straightforward to 
prove [3^ 1 that the eigenvalues of a circulant matrix 

/ c ci ■ • ■ c/c-i \ 



C 



c ci 

Cfc-i c 



C-2 



Cfe-2 
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are given by cq + c\Tj + • • • + Cfe_ir| _1 , where rj = 
exp(2-7rij / TV) is an TVth root of unity. Thus, the eigen- 
values /3j{n) of J(n) are given by 

&(n) = Mn((l- £ ) + |(r,+rf- 1 )) 
= Vn<Pj(e,N), 

where 

i (e,JV) = (l-e)+ E cos^V (7) 

It is important to notice that ^ (e, iV) does not depend 
on the iteration n: the time dependence has been de- 
coupled (factorized) into \x n . The Lyapunov exponents 
are then given by 

t 

A, = lim InTT \Pi{n)\ l,t 

n=l 

= Mm In M^iCe, JV)| Ia*"I 1A J 
1 * 

= ln|&(e,JV)l + lim - Vln|/i„|. 

t— »oo t * — ' 

n— 1 

Thus by defining Ao to be the Lyapunov exponent of a 
typical orbit of a single, uncoupled, local map, starting 
at a; : Ao = limt_ >00 (l/f) J2n=i m \^n\, one obtains the 
following expression for the Lyapunov exponents of a 
homogeneous evolution: 

Aj = Ao + In \4>i(e, N)\. (8) 



derivative of the local map at any point is constant. 
Examining further this similarity, if one takes the fully 
chaotic logistic map (4 x (1 — a;)) as the local map for the 
diffusive CML, the LS for the homogeneous evolution is 

A, : = ln2 + In \<pi(e,N)\. (10) 

In fact, any one-dimensional map whose Lyapunov ex- 
ponent is Ao = In 2 gives rise to the LS ( |Io| ) under homo- 
geneous evolution. The LS ([!(]) corresponds exactly to 
the LS of a lattice of coupled Bernoulli shifts and thus 
the results described below for the rescaling of the sub- 
system LS are valid for the case of a lattice of coupled 
Bernoulli shifts. 

Now let us perform the LS analysis for a sub-system of 
the original CML. Thus instead of taking all the sites i — 
1, . . . , N we take N s sites starting at any position j. The 
choice of j is not important since we are dealing with 
periodic boundary conditions and because the state is 
homogeneous; from now on we choose j = 1. Thus, we 
take a principal sub-matrix J' of size N s x N s from the 
whole Jacobian J. In matrix terms J' = 7r( J) where 7r 
is the following projection 

tt( J) = IT • J • IL. (11) 

with the left (11;) and right (II r ) projection matrices 
defined as 

IL = (l\z) 




Note that the Lyapunov exponents defined by (JSJ) are 
not arranged in decreasing order. Re-indexing them in 
decreasing order they become 

!A + In |0 fe (e, N)\ k even 
(9) 
A + In \<j)k=i{e,N)\ k odd 

where k = 1 to N. It is clear that Xk = Xk+i when k is 
even, so most of the exponents occur in degenerate pairs, 
apart from the largest, and, if N is even, the smallest. 
The linear stability of a homogeneous orbit is then char- 
acterized by the Lyapunov exponent Ao of a single site 
in the uncoupled case (e = 0). In particular, if the local 
map is not chaotic then the homogeneous evolution is 
not chaotic either since Afc < Ao (|</>fc(e, iV) | < 1 for all 
k). 

It is interesting to notice that the same shape for the LS 
of a homogeneous CML (cf. (pf) ) is obtained for a lattice 
of coupled Bernoulli shifts |33j for any orbit. There is 
however an important difference: while in the CML the 
LS dependence on the actual orbit was decoupled thanks 
to the homogeneity, in the case of coupled Bernoulli 
shifts, the LS is decoupled from the orbit because the 



where, from now on, I is the iV s x N s identity matrix 
and Z is the N s x (N — N s ) null matrix. Therefore, in 
order to compute the Lyapunov exponents for the trun- 
cated system one has to compute the following product 
of projected matrices 

P'{n) = (UiJ(n)U r )---{UiJ(2)U r )(IliJ{l)Ilr) 

(12) 

= 11* J{n)TL c J{n -!)■■■ J(2)U c J(l)U r , 

where 





z 








Multiplying equation ( |l2|) from the left by the identity 
matrix obtained by II; • II r = / yields 

p\n) = n z n r iLJ(n)n c j(n-i)--- j(2)n c j(i)n r 
= n ; n c j(n)u c j(n - 1) ■ ■ • j(2)n c j(i)n r 

= (P(n)) , 

(13) 
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where we 
define the new product P(n) = K(n) ■ ■ ■ K(2)K(1) of 
the projected matrices K(i) = H c J(i). 

Using the above description we obtain the sub-system 
LS for the homogeneous evolution. The projected Jaco- 
bian for the homogeneous evolution at time n is 

J'(n) = 7r(J(n)) = /i„7r(M) =fi n -M' (14) 

where M' = ir(M) is the N s x N s constant matrix 

/1-e e/2 •■• \ 
e/2 1-e e/2 ■■■ 
M / = e/2 1 - £ • • • 

V • • • e/2 l-e J 

if 7V S < and Mjy = M if iV s = iV. From now on 
we only use the notation M' when N s < N. It is im- 
portant to notice that by taking a sub-Jacobian matrix 
the periodicity of the boundary conditions is lost. The 
dynamics of the sub-system at the boundaries could be 
thought as being coupled to some external noise com- 
ing from the adjacent sites. Thus, in contrast to M, the 
matrix M' is not circulant, however its eigenvalues are 
well known to be |34| 

0;( £ ,jV s ) = (l-e)+£cos(-^). (15) 

(where j = 1 to N s ); and so the eigenvalues fi'Ari) of 
J'{n) are /3'-(n) = ix n 4>'As,N s ). The sub-system LS is 
given by 

A; = A +ln|^.( e ,7V s )|. (16) 

One can immediately infer from this that the Lya- 
punov exponents for the homogeneous evolution are in- 
terleaved for two consecutive sub-system sizes. More 
precisely, suppose that we take two sub-systems, one of 
size N 8 and the other of size N s + 1. It is then trivial 
to see that their respective Lyapunov exponents A^(iV s ) 
and A-(iV s + 1) satisfy: 

K(N S + 1) < \i(N s ) < K+i( N s + 1) VI < i < N s , 

(17) 

see figure |l|.a. Interleaving of the sub-system LS with 
respect to the whole LS Aj(iV) also occurs: 

\i{N) < X^Ns) < \ l+N - Ns (N) Vl<i<N s . 



the interleaving property of the Lyapunov exponents for 
the homogeneous case is a straightforward consequence 
of the decoupling of the time dependence of the Jaco- 
bian matrix leaving us with the constant matrices M 
and M 1 . In a typical non-homogeneous evolution the 
time dependence of the Jacobian cannot be factorized 
and an equivalent constant matrix for the Jacobian does 
not exist. Therefore, Cauchy's interlace theorem cannot 
be applied in this general case and there is no reason a 
priori for the interleaving property to hold for a generic 
extended dynamical system. It is true that, at any par- 
ticular time, there is interleaving between the eigenval- 
ues of the whole Jacobian and those of a sub-system. 
However, when computing the LS, one has to compute 
the product of the Jacobian matrices while for the sub- 
system LS one uses the product of the sub-Jacobian 
matrices and therefore the interleaving of the matrix 
product is no longer assured. The only way, a priori, for 
the interleaving to work would be to take the product of 
the whole Jacobians first and only then extract the sub- 
Jacobian. The problem with this procedure is that one 
has to rely again on re-orthonormalization procedures 
involving the original matrix size N, making the task 
impossible for large N. Nevertheless, as we shall see in 
the following section, the interleaving of the sub-system 
LS does hold to a great extent in the thermodynamic 
limit. 




This interleaving of the eigenvalues is a consequence of 
Cauchy's interlace theorem [19| that gives bounds on the 
eigenvalues of a principal sub-matrix given the eigenval- 
ues of the original matrix. It is important to notice that 
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FIG. 1. Lyapunov spectrum for a homogeneous evolu- 
tion in a diffusive CML: a) interleaving for sub-system sizes 
N s = 1, . . . , 20 (N = 20); b) rescaled sub-system LS, the cir- 
cles represent the whole LS (iV = 30), the thin dashed lines 
represent the functions A dd and A cvcn passing through the 
eigenvalues for even and odd indexes respectively, while the 
thick lines represent the rescaled LS with N s — 10 using the 
conventional rescaling r' = N/N s (thic k dashed line) and the 



new rescaling obtained in section 
(thick solid line) 



II A 



r = {N + l)/(N s + l) 



Another important point to note from equations ( [lq ) 
and ( |l6| ) is that the LS of the sub-systems all have the 
same shape. The best way to see this is to rescale the 
indices of the Lyapunov exponents so that they lie in the 
range [0,1]: so instead of plotting A against j we plot it 
against j/(N s + 1). Equations (15) and ( |l6| ) then show 
that the points always lie on the graph of the function 



X(z) = A + In [(1 — e) + e cos(7rz)] 



(18) 



irrespective of the value of N s . This observation sug- 
gests another way of looking at the interleaving prop- 
erty. For a given N s , the z values of the sub-system LS 
are equally spaced in the interval [0,1]; if we increase N s 
by 1 the new z values interleave with the old. Since A is 
a monotone function the fact that the z values interleave 
means that the X(z) values interleave also. It is worth- 
while mentioning that we are considering the simple case 
1 — 2e > so the absolute value inside the logarithm 
in equation ( |l6| ) can be omitted. For 1 — 2s < the 
eigenvalues need further re-indexing in order to main- 
tain their decreasing order and a similar construction as 
bellow is possible. 

To compare the sub-system LS with that of the full sys- 
tem we should similarly rescale the indices for the lat- 
ter, so now we plot the full system Lyapunov exponents 
against j/(N + 1) instead of j. The points of this spec- 
trum do not lie on the graph of A(z); however, equation 
@ shows that they do lie on the graphs of the functions 

\even(z) = A (z (1 + 1/N)) 

(for exponents with even indices) and 



A o dd(z) = A(z(l + l/A0-l/A0 

(for exponents with odd indices), where the function 
A(z) is given by equation (Jisj). Since z(l+l/N)-l/N < 
z < z(l + 1/N) (0 < z < 1) and X(z) is a decreas- 
ing function we see that X cvcn (z) < X(z) < X Q dd(z)- 
Thus Aovon and A dd are bounding curves for A (see 
thin dashed lines in figure |.b) and converge to it as 
N — > oo; the differences between A and the other curves 
are 0(1/N). 

The similarities between the shapes of the Lyapunov 
spectra of the sub-systems and of the whole system 
mean we can use the sub-system LS to estimate the 
whole LS: to do this we rescale the indices of the sub- 
system exponents, plotting Xj against rj where r is a 
factor chosen so that the rescaled sub-system LS lies as 
close as possible to the plot of the full system LS. The 
above discussion shows that if we choose 



N+l 



1 



(19) 



then the rescaled sub-system LS differs from the full 
system LS by an amount of 0(1/N). 

The scaling given by (|l9|) differs from that used conven- 
tionally, which is performed by scaling by 

r> = — 

N s 

see p|Jl^ |l^ , ^5|] . It is clear however that using r' will 
give results that differ from those using r by terms of 
0(1/N S ), and since this is larger than 0(1/N) the er- 
rors in the exponents will also be 0(1/N S ). This sug- 
gests that scaling (|l^) should give more accurate results 
than the conventional scaling; this is certainly true in 
the homogeneous case. As an example figure 0.b shows 
the original LS for a homogeneous CML with N = 30 
(circles) along with the rescaled LS with N s = 10 using 
the conventional rescaling r' (dashed line) and the new 
rescaling r obtained above (solid line). It is clear that 
the new rescaling gives a much better approximation to 
the original LS. 



B. Interleaving and rescaling for coupled logistic 
maps 

As mentioned in the previous section, the interleaving 
property for the homogeneous evolution relies on the 
fact that the Jacobian matrices can be factorized into 
a time dependent scalar and a time independent matrix 
(see equations @ and (fl4|)). For a non-homogeneous 
evolution the Jacobians cannot be factorized in such a 
way and thus a priori one does not expect interleav- 
ing to occur. Surprisingly enough the numerical evi- 
dence points towards interleaving of the sub-system LS 
for almost every Lyapunov exponent in the fully devel- 
oped chaotic regime. In this section we shall present 
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such evidence for a logistic CML, and discuss why such 
behaviour might be expected to occur. More general 
systems will be considered in the following section. 

We thus consider the diffusive CML (|j) with the fully 
chaotic logistic map f(x) = Ax (1 — x) and compute its 
LS for several values of the coupling parameter e. As 
with all numerical work in this paper we employ a fast 
HQR algorithm for the computation of Lyapunov expo- 
nents ||. We then calculate the sub-system LS using 
principal sub-matrices J' of size N s — 1, ... ,30 of the 
Jacobian. In doing so one is not taking into account the 
dynamics of the neighbouring sites next to the boundary 
and their effects are considered as noise. Thus, the algo- 
rithm consists in computing the LS of the sub-Jacobian 
J' by truncating the actual Jacobian J at each time step 
and then applying the HQR algorithm. The results are 
shown in figure where we plot the sub-system LS for 
increasing sub-system size (N s — 1, . . . , 30) for 3 differ- 
ent values of the coupling parameter. In the figure, the 
filled circles represent the Lyapunov exponents that do 
not fulfil the interleaving condition. Strikingly, the LS 
corresponding to e — 0.05 and e = 0.45 (figures a and b) 
are very well interleaved, with the exception of a couple 
of points. On the other hand, for e = 0.95 (figure c) the 
LS is not that well interleaved for the smallest Lyapunov 
exponents, although for the large ones the interleaving 
is as good as for the previous two figures. The reason 
for this failure for the smallest Lyapunov exponents is 
that in the limit e — ► 1 the lattice decouples into two 
independent sub-lattices: one for odd i and the other 
for even i. Thus, when successively increasing the sub- 
system size, one is including in turn contributions from 
the even and the odd sub-lattice. This is reflected in 
a variation in the smallest Lyapunov exponents every 
time we increase the sub-system size by one, hence the 
bi-periodic nature of the interleaving failure. In fact, 
by removing the sub-system LS for odd sizes one ends 
up with almost perfect interleaving. The exact reasons 
and conditions for the interleaving of the sub-system LS 
to happen are not yet understood, however we believe 
that they are connected with the convergence of the sub- 
system LS to the full system LS — a convergence which 
may be expected in the thermodynamic limit, see below. 
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FIG. 2. Sub-system Lyapunov spectra for the fully chaotic 
coupled logistic lattice N = 100 for sub-system sizes 1 to 30 
(left to right) for a) e = 0.05, b) e = 0.45 and c) e = 0.95. 
The filled circles represent those Lyapunov exponents which 
fail to interleave. 
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FIG. 3. Comparison of the whole Lyapunov spectrum 
(solid line) and the rescaled sub-system Lyapunov spec- 
trum using the new rescaling r (circles) and the conventional 
rescaling r' (crosses) in the fully chaotic logistic lattice with 
N = 100 for several sub-system sizes (N s — 15, ... , 25). a) 
e = 0.05 and b) e = 0.45. 

As mentioned in the introduction, it has been observed 
for some time that under appropriate rescaling the sub- 
system LS approximates the whole LS. The usual argu- 
ment for this rescaling behaviour makes use of the ther- 
modynamic limit. In the previous section, while study- 
ing the interleaving of sub-system LS for the homoge- 
neous case, a new rescaling was suggested (see equa- 
tion (fl9|)). Let us test this for the case of the fully 
chaotic coupled logistic lattice. In figure || we com- 
pare, for e = 0.05 and e = 0.45, the rescaled sub-system 
LS using the new rescaling r = (N + l)/(N s + 1) © 
(circles) and the conventional one r' — N/N s (crosses) 
to the whole LS (lines) for different sub-system sizes 
(N s = 15, . . . , 25). As is clear from the figures, the new 
rescaling r gives a much better fit to the original LS 
than the conventional rescaling. 
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FIG. 4. Rescaled Lyapunov spectrum (circles) for the 
coupled logistic lattice with e = 0.95 for sub-system sizes 
N s = 15, ... , 30. The solid line represents the whole Lya- 
punov spectrum N s — N = 100. 



Let us explore the idea of rescaling the sub-system LS 
in the thermodynamic limit a bit further. The corre- 
spondence between the rescaled LS and the whole LS in 
figure H is astonishingly good. The rescaled spectra lie 
almost perfectly on top of a decreasing curve, therefore, 
as with the homogeneous case discussed above, it is not 
surprising that they are interleaved. In general, if the 
rescaled Lyapunov spectra of the sub-systems converge 
sufficiently quickly to the whole system LS we expect 
to have good interleaving of the sub-system Lyapunov 
spectra. On the other hand, if the rescaled sub-system 
LS do not approximate the whole system LS well, it is 
not clear that interleaving will occur. To illustrate this 
we present the rescaled LS using the new rescaling r for 
e = 0.95 in figure |^. In this case, the rescaled LS do 
not give such a good approximation to the whole LS (in 
particular for the smallest Lyapunov exponents) as seen 
in the other cases (e = 0.05 and s = 0.45). As explained 
above, this is due to the decoupling of the whole lattice 
into two sub-lattices when e — ► 1. Therefore it appears 
that the non-interleaving of the smallest Lyapunov ex- 
ponents in figure ||.c is related to the lack of convergence 
of the sub-system LS. In general we suppose that fail- 
ure to interleave is an indication that the sub-system 
LS have not converged. Clearly however, the presence 
of interleaving is not a sure indication that convergence 
has occurred; this is illustrated by the two-dimensional 
logistic lattice discussed below. 

We believe that the key point in understanding the in- 
terleaving behaviour is that although in computing the 
sub-system LS one is using the product of projected ma- 
trices ([l3]), one does not modify the original dynamics 
in any way. Recall that similar matrices share eigen- 
values. Thus a feasible explanation for the occurrence 
of interleaving is to hypothesize that the product of the 
projected matrices P'(n) is a projection of a N x TV ma- 
trix Q(oo) which is similar to the limit as n — > 00 of 
the original product P(n) of the whole Jacobians. In 
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other words, we conjecture that there exists an invert- 
ible N x N matrix S such that 



Q(oo) = lim S^P^S, 



(20) 



where the product of the projected matrices P'(n) in 
the limit is obtained by projecting Q(oo): 

P'(oo) =tt(C(oo)). 

Informally, this is saying that in some sense in equation 
( ^3| ) the projection matrices commute on average with 
the Jacobians in the n — > oo limit. We believe that it 
might be possible to make this statement rigorous by an 
appropriate generalization of the multiplicative ergodic 
theorem. 
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FIG. 5. Interleaving of the sub-system LS, N s = 1, . . . , 30, 
for the fully chaotic logistic lattice with e = 0.45 using 
the more general projection matrices (|^) to extract the 
sub-system Jacobians: a) II i and b) II2. 



One might then also ask what is so special about the 
projection II C . Is it possible for interleaving to occur 
for more general projections? The following two exam- 
ples suggest that this is indeed the case. Consider the 
following projection matrices 



n, = 



n 2 = 



/ , 
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y Z tr 




1 n 2 
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(21) 



where Z is the N s x(N~N s ) null matrix and the N s x N s 
matrices LT^ and W 2 are 



m = 
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FIG. 6. Rescaled sub-system LS corresponding to figure pi 
using the projection matrices IIi (crosses) and IJ2 (circles J. 
The continuous line corresponds to the original LS computed 
with the whole Jacobian. 



where the are random numbers chosen from the in- 
terval [0, 1] with equal probability. Note that we are 
still using the term projection matrices for ill and II2 
which in a strict sense is not correct, since they do not 
satisfy ELj = Il| (j = 1, 2). We use this terminology to 
stress the fact that they completely remove some of the 
entries of the original Jacobian. Thus, instead of taking 
the projection matrix LI C let us take 111 and II2. For 
the projection II2 we randomise its entries every time- 
step; similar results were obtained by randomising only 
at the beginning and keeping the same projection ma- 
trix thereafter. In figure @ we depict the non-rescaled 
sub-system LS using both projection matrices for the 
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fully chaotic logistic lattice. The figure strongly sug- 
gests that in these cases interleaving still occurs. It 
thus appears that the choice of projection matrix is not 
a crucial ingredient for interleaving. Nonetheless it is 
important to say that we do not expect interleaving to 
hold if one uses a series of projection matrices such that 
when computing the LS one does not get convergence. 
In the above examples, III and II2, we do have the re- 
quired convergence. For the II2 case, the convergence 
of the LS of their product is a well known fact ]3q]37[| . 

On the other hand, when we turn to rescaling we find 
that although for III and II2 we still get convergence 
of the rescaled sub-system LS to a definite limit, this 
limit is not the original LS for the full system (figure ^) . 
The reason for this discrepancy is easy to understand 
since the new projections IT and II2 combine the entries 
of the projected Jacobians and thus one expects the 
eigenvalues to change. 



Estimation of quantities derived from the 
Lyapunov spectrum 



As illustrated in the previous section, the LS can be 
well approximated by the rescaled sub-system LS in the 
thermodynamic limit. We now use the new rescaling in 
order to approximate the original LS by extrapolating 
from the sub-system LS. We estimate the largest Lya- 
punov exponent, Lyapunov dimension and KS entropy 
and we compare our method to the results obtained with 
the whole LS and with the conventional rescaling. 

The first method to approximate quantities derived 
from the LS in the thermodynamic limit is by defining 
intensive quantities from the extensive ones by simply 
using the corresponding densities [38|Jl^Jl^ , p5| ,p[ . Let us 
define the densities of ([!]) and (|J): 



Pd(N s ) 
Ph(N s ) 



El 

h 



(22) 



corresponding to the Lyapunov dimension density and 
the KS entropy density respectively. In the thermo- 
dynamic limit these densities are intensive quantities 
(i.e. they do not depend on the sub-system size taken). 
One then estimates their extensive counterpart when 
N s — > N by multiplying the densities ( p2] ) by N. To es- 
timate the largest Lyapunov exponent for the whole sys- 
tem we directly take the value of the largest Lyapunov 
exponent of the sub-system (the Lyapunov exponents 
are not extensive quantities). It is worth mentioning 
that in order to use these intensive densities to estimate 
extensive ones we are supposing the size N of the orig- 
inal system to be known. 




90 



D L (N)- 
rescaling r — e— 
rescaling r' — 
N x /> (/ (/Y s ) •> 



<-x-x-x-x-x-x-> 



b) 



20 



7 - 



h(N) 

rescaling r — e— 
rescaling r' — 
Nxp h {N s ) -o- 




j-o-o o -0-^0- o o -0-0- o o -0-0- ©-00 -o-o- o 



c) 



5 10 15 20 25 30 35 40 

N s 

FIG. 7. Estimation of a) the largest Lyapunov exponent, 
b) the Lyapunov dimension and c) the KS entropy as a func- 
tion of the sub-system size N a in the coupled logistic lattice 
with TV = 100 and e = 0.45. The estimates obtained by 
using a) the largest Lyapunov exponent of the sub-system 
and b-c) the associated densities from the sub-system are 
presented with diamonds, and the estimate obtained from 
the piece-wise linear fit to the rescaled LS is presented with 
crosses for the conventional rescaling and circles for the pro- 
posed new one. The values obtained with the whole LS are 
represented by the horizontal line. 
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The second method, which we believe is more accurate, 
consists of taking the sub-system LS, rescaling it, ex- 
trapolating a curve through it to obtain an approxi- 
mation to the whole LS and only then computing the 
desired quantities. There are several ways to extrap- 
olate the whole LS from the sub-system LS; here we 
have chosen a piece-wise linear approximation for sim- 
plicity. One could use more accurate methods such as 
cubic splines but the aim here is to compare both kinds 
of rescaling and thus a piece-wise linear fit is the most 
straightforward approach. Therefore, take the rescaled 
LS Xi(N s ), obtained with either rescaling for a sub- 
system of size N s , and consider the polygon V through 
all the points (i, Xi(N s )). To estimate a Lyapunov expo- 
nent of the whole LS lying between Ai (N s ) and Xn„ (N s ) 
one simply uses the fit given by the polygon V . For Lya- 
punov exponents lying to the left (right) of the polygon 
use linear extrapolation from the first (last) two points 
of the rescaled LS. Here again one could use more so- 
phisticated extrapolation methods but for simplicity we 
restrict ourselves to the linear one. Once the whole LS 
is estimated using the above method, or a more com- 
plicated one, quantities such as the largest Lyapunov 
exponent Ai(iV), the Lyapunov dimension Dl and the 
KS entropy h are easily extracted. 

In figure |^ we compare the estimates of a) the largest 
Lyapunov exponent X\(N), b) the Lyapunov dimension 
Dl and c) the KS entropy h obtained from the in- 
tensive densities (diamonds) and the piece-wise linear 
fitting for both rescalings (conventional rescaling with 
crosses and the proposed new one with circles) as the 
sub-system size increases for the coupled logistic lat- 
tice. The actual values of these quantities calculated 
with the whole LS correspond to the horizontal lines. 
For the largest Lyapunov exponent, figure Q.a, we notice 
that the estimates are almost identical for both rescal- 
ings (crosses and circles). This is due to the fact that 
both rescalings tend to coincide for small i (see figure 
0.b). The estimate of the largest Lyapunov exponent by 
just taking the largest Lyapunov exponent of the sub- 
system (diamonds) shows a slower convergence than the 
linear fit methods. For the Lyapunov dimension, fig- 
ure ^.b, the method with the slowest convergence cor- 
responds to the conventional rescaling (crosses), while 
the approximations derived from densities (diamonds) 
and from a linear fit with the new rescaling (circles) are 
quite good (note that the new rescaling method does 
better than the approach using densities). Finally, for 
the KS entropy, figure [?].c, the estimates using the den- 
sity (diamonds) and the conventional rescaling (crosses) 
have similar convergence rates while the new rescaling 
method (circles) does considerably better. The evidence 
given by this set of plots tends to indicate that the new 
rescaling method gives better convergence to the quan- 
tities derived from the sub-system LS. 



III. MORE GENERAL EXTENDED 
DYNAMICAL SYSTEMS 



So far we have only considered interleaving and rescaling 
in systems in one spatial dimension with nearest neigh- 
bour coupling, corresponding to tridiagonal Jacobians. 
In this section we turn to more general kinds of extended 
dynamical systems by allowing a larger coupling range 
(e.g. chaotic neural networks) and by taking a different 
topology for the lattice (e.g. lattice with two spatial di- 
mensions). The results presented in this section suggest 
that the interleaving and rescaling properties observed 
for the simpler one-dimensional CML persist for more 
general extended dynamical systems. 



A. Chaotic neural networks 



We now consider a chaotic neural network jl2| of the 
form 



x? +1 = tanh 




(23) 



where g is a real number called the gain parameter, 
k represents the connectivity (essentially playing the 
same role as the range of the coupling in a CML) 
and the weight matrix Cy has entries chosen randomly 
from [—1,1] with a uniform probability distribution for 
(i — j) (mod N) < k and CV,- = otherwise. 

Both the CNN and CML dynamics work in two stages — 
nonlinearity and coupling — but their order is inverted. 
The CML dynamics applies the nonlinear mapping / 
first and then the coupling, while the CNN first applies 
the coupling via a linear weighted combination of neigh- 
bouring sites, and then a nonlinear map (the sigmoid). 
This inversion is reflected in the Jacobian matrix of the 
transformation: while each entry of the CML Jacobian 
(||) depends on a single site, each entry of the CNN 
Jacobian depends on a neighbourhood of sites: 



Jij(n) 



gCi. 



cosh 



i+k 



1= 



k Cil x ? 



(24) 



The CNN Jacobian (|24|) inherits the zeros of the cou- 
pling matrix Cy, i.e. Jij(n) = if (i — j) (modiV) > k. 
Another difference between the CNN that we will con- 
sider and the diffusive CML discussed before is that 
the CNN involves coupling with a larger neighbourhood 
than just the left and right nearest neighbours. 

Let us now analyze the interleaving and rescaling for a 
CNN with a large k. In figure || we show the interleaving 
and rescaling with k = 10 and g = 2. As we can see, 
the interleaving is quite good with the exception of a 
few small Lyapunov exponents. In figure |^.b we plot 
the rescaled sub-system LS for several sub-system sizes 
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using both rcscalings (circles: new rescaling and crosses: 
conventional rescaling) along with the whole LS (solid 
line). Clearly the new rescaling gives a better estimate 
of the whole LS. Similar results were obtained for other 
values of the parameters k and g. 
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FIG. 8. a) Interleaving of the sub-system LS in the chaotic 
neural network (E3I) with k = 10 and g = 2. b) Comparison 
between the conventional rescaling of the sub-system Lya- 
punov s pect rum (crosses) and the new rescaling obtained in 
section [I A (circles), the whole LS is depicted by the solid 
line. 



B. Two-dimensional logistic lattice 



The interleaving and rescaling proper ties of the sub- 
system LS were obtained in section II A for a one- 
dimensional array of coupled maps. Here we put to the 
test the interleaving and rescaling for a two-dimensional 
CML. Let us take a two-dimensional square lattice of 
size L x L. The local dynamics x™ at each node 
and any time n is governed by the fully chaotic logistic 
map 



As in the one-dimensional CML the local dynamics is 
applied first 

y%{x) =/K), 
and then the coupling dynamics 

4 +1 = (i- £ )^ + ^., 

where . is the average of the j/™- in the neighbourhood 
Mij of site (i, j). The neighbourhood A/y is taken to be 
the eight adjacent sites to with periodic boundary 
conditions. 

The Jacobian J{n) at time n for this two-dimensional 
lattice is defined through its elements: 

where the indices <Jk and o~i refer to the position in the 
actual lattice of the chosen fcth and Zth state variables 
of the system. If one just wants to compute eigenval- 
ues of the whole Jacobian, the order in which the state 
variables are taken is not relevant. However we are in- 
terested in extracting sub- Jacobian matrices from the 
whole system and thus the ordering choice of the state 
variables does matter. There are L 2 \ different ways to 
choose the ordering, but the simplest way consists of 
taking the site (1,1) as the first state variable and then 
proceeding horizontally to the right until the end of the 
lattice is reached and then proceeding to the bottom of 
the lattice by rows: 
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fij(x) = f(x) = Ax (I - x). 



that is <Tk = (k — [k/L\,[k/L\) where [z\ denotes 
the largest integer smaller than or equal to z. From 
now on this kind of ordering will be called horizontal 
wraparound. There is obviously a vertical counterpart 
where the order is taken by columns. The problem with 
this type of ordering is that it does not build up the 
Jacobian in a natural way. The propagation of a per- 
turbation typically grows equally in both of the two 
spatial dimensions (in particular for our choice of cou- 
pling since all the neighbours contribute with the same 
weight e/8). In contrast, with horizontal or vertical 
wraparound one has to wait until a complete wrap is 
taken to fall again near the perturbed area. A more 
natural approach might thus be to attempt to mimic 
the spatial growth of perturbations by taking an order- 
ing that fills up a two-dimensional area from the centre 
outwards. For that purpose, we use the following order- 
ing technique: 
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We call this square wraparound. 
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FIG. 9. Sub-system Lyapunov spec- 

tra for the two-dimensional 20 x 20 coupled logistic lattice 
for sub-system sizes 1 to 40 (left to right) for e — 0.45 and 
for the two wraparound methods for building up the Jaco- 
bian: a) square wraparound and b) horizontal wraparound. 
The filled circles represent the Lyapunov exponents where 
interleaving fails. 



in the computation of the exponents (and in particu- 
lar poor convergence). Therefore, we shall consider a 
Lyapunov exponent to be interleaved if it falls in the 
interval defined by the inequality (|i~7| ) with an error 6: 

X t (N s + 1) - SA < Xi{N s ) < X l+1 (N S + 1) + 5A, (25) 

where A = Xi + i(N s + 1) — Xi(N s + 1). From now on we 
redefine 5 such that the errors are given in percentages. 
Using such a definition, if one allows a small error of 
2.5% — S = 0.025 in j25| ) — for the Lyapunov exponents 
in figure ^, one obtains perfect interleaving for the whole 
spectrum. 
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FIG. 10. Rescaled sub-system LS for the two-dimensional 
coupled logistic lattice (same parameters as in figure |^) using 
a) square and b) horizontal wraparound methods. 



In figure |9| we show the non-rescaled sub-system LS for 
the two wraparound methods, a) square and b) horizon- 
tal, and we plot with solid circles the Lyapunov expo- 
nents that fail to interleave. Observe that interleaving 
failure occurs for only a very few Lyapunov exponents. 
After a careful examination of these Lyapunov expo- 
nents one notices that they are very close to interleav- 
ing, suggesting that the failure is due to numerical error 



The interleaving seen in figure |9j suggests that the or- 
dering choice for the Jacobian entries does not play an 
important role in this phenomenon. However, as can 
be seen in figure [Io[ where we plot the rescaled LS for 
both wraparound methods along with the whole LS, the 
choice of ordering method is crucial in obtaining good 
rescaling behaviour. Square wraparound (figure |i"(i|.a) 
yields immediate convergence towards the whole LS: 
even for a very small sub-system size the rescaled LS 
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is almost exactly superimposed on top of the whole LS. 
On the other hand, horizontal wraparound (figure p"ojb) 
gives a rescaled sub-system LS that seems to converge 
to a different curve for sub-system sizes N s = 1, . . . , 20 
(aligned points in the lower part of the curve for the 
first Lyapunov exponents). For sub-system sizes larger 
than 20 the rescaled LS starts a new convergence to- 
wards something closer to the whole LS. The explana- 
tion for this phenomenon is quite simple. The Jacobian 
for the horizontal wraparound consists of a main diago- 
nal of non-zero elements coming from the neighbours in 
the same row of the square lattice, however, the neigh- 
bours in the row above and below give rise to two sub- 
diagonals of non-zero elements. The sub-diagonals start 
when a whole wraparound has been completed, that is 
when N s — L where L is the side length of the square 
lattice. Thus for sub-system sizes N s < L the sub- 
Jacobian only extracts the main diagonal elements and 
does not capture the two sub-diagonals with vital infor- 
mation about the neighbouring sites in the rows above 
and below. When N s > L the sub- Jacobian starts cap- 
turing these forgotten neighbours and the rescaled sub- 
system LS now begins to converge to the desired LS. 
For the example in figure [l0].b this behaviour starts at 
N s = L = 20. This effect of horizontal wraparound is 
reflected when one tries to extract information from the 
sub-system LS. As an example, we depict in figure [ll] an 
estimate of the largest Lyapunov exponent by extrapo- 
lating the whole LS from its rescaled version as the sub- 
system size increases. The results are depicted with cir- 
cles for square wraparound and with crosses for horizon- 
tal wraparound. The vertical solid line corresponds to 
the largest Lyapunov exponent from the whole LS. The 
estimate using horizontal wraparound seems to converge 
to a much smaller value than the desired one for sub- 
system sizes N s < L — 20. When the sub-system size is 
increased further, horizontal wraparound performs bet- 
ter but still lacks the desired convergence. On the other 
hand, the square wraparound converges rapidly in a 
smooth way: this is because it was designed to build up 
the Jacobian entries in a more natural way. Therefore, 
although the interleaving for both wraparound methods 
is very good it is considerably more reliable to use the 
square wraparound for rescaling purposes of the sub- 
system LS. 
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FIG. 11. Estimation of the largest Lyapunov exponent 
as a function of the sub-system size N s in a two-dimensional 
logistic lattice of size 20 x 20 and with e = 0.45 using a linear 
fit for the rescaled sub-system LS. The circles correspond 
to building up the Jacobian by square wraparound whilst 
the crosses correspond to horizontal wraparound. The value 
of the largest Lyapunov exponent for the whole lattice is 
represented by the horizontal solid line. 



C. Host-parasitoid system 



We now consider a more general type of two dimen- 
sional lattice, namely the Host-Parasitoid lattice model 
[ p9| [II]] . For this system the local dynamics is no longer 
one-dimensional but two-dimensional: hosts and para- 
sites. The model evolves again in two phases. First 
there is at each site a local dynamics given by 



CH Z ( 



1 



-a PI 



(26) 



where H n and P n are respectively the population size 
of hosts and parasitoids at time n, a is the per capita 
parasitoid attack rate, b is the host reproductive rate 
and c is the conversion efficiency of parasitised hosts into 
female parasitoids in the next generation. The second 
phase involves dispersal into a neighbourhood A/y of site 
i.e. a fraction fi^ of hosts and \x v of parasitoids 
disperse equally into the eight neighbouring sites: 



(27) 



where "Hj^ i . and Vj^ i . are, respectively, the average of the 



hosts and the parasitoids (after local dynamics (26)) in 
the neighbourhood A/y of site We take a square 

lattice € [1,£] 2 and periodic boundary conditions. 
The total size of the system is then TV = 2L 2 . Let us 
build up the whole Jacobian with host-parasite blocks 
of size 2x2: 
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where the host-parasite blocks J£* are given by 
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FIG. 12. Interleaving of the sub-system LS for the 
host-parasite system in a two-dimensional lattice of size 
20 x 20. The Jacobian was built using a)-b) square 
wraparound and c)-d) horizontal wraparound. Figures b) 
and d) correspond, respectively, to amplifications of figures 
a) and c) for the top half of the spectrum. 



The indices <j\ l 2 refer to the actual position in the 
two-dimensional lattice of a particular local population. 
As for the two-dimensional CML the ordering choice 
of the Jacobian entries plays an important role for the 
rescaling. 

Given a reasonable lattice size (L > 15) and depending 
on the dispersal parameters fih and fj, p the evolution of 
model ( p7| ) is spatio-temporally chaotic p2j. Here we 
choose L = 20, a = 1, b = 2, c = 1, fih = 0.2 and 
fi p = 0.6. The full system is thus N — 2L 2 — 800 di- 
mensional. We start the system with random initial con- 
ditions and discard a transient of 10 5 iterations before 
computing the sub-system LS. In figure |l2| we depict 
the interleaving of the sub-system LS for sub-system 
sizes N s = 1, . . . , 40 where we allow a 5% error in the 
interleaving — S = 0.05 in (^5|). Figures [l^.a and p^[b 
correspond to square wraparound whilst figures p^.c and 
[l2|.d correspond to horizontal wraparound. As the figure 
shows, interleaving is quite good even for the upper re- 
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gion (see amplifications in figures [lj.b and |l2].d) where 
the density of Lyapunov exponents is very high and the 
intervals for interleaving are small and thus the mar- 
gin for error in the inequality ( p5| ) is reduced. Square 
wraparound does better for large Lyapunov exponents 
(figure |l2|.b) whilst horizontal wraparound does better 
for small ones. However, overall both methods have ap- 
proximately similar performance. 
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does square wraparound. The reason is again that for 
the horizontal wraparound one has to wait until a com- 
plete wrap is finished until falling again into the neigh- 
bouring region. In this horizontal wrap of the 
Jacobian is achieved when N s = 2L = 40. Moreover, 
for N s — 2L = 40 one is only including partial deriva- 
tives of hosts with respect to hosts and parasitoids. In 
order to include dependences of parasitoids with respect 
to hosts and parasitoids one should take a further wrap 
of the Jacobian, i.e. N s = AL — 80. 

Therefore, the horizontal wraparound technique for sub- 
system sizes N s < 40 does not pick up the dynamics of 
the neighbours situated in adjacent rows. This problem 
for horizontal wraparound becomes worse as the dimen- 
sion of the local dynamics is increased. A partial solu- 
tion to this problem is to build up the Jacobian by using 
just one of the local variables of the system. Particu- 
larly in the host-parasite system where the parasitoid 
dynamics is slaved to the host dynamics, one should 
be able to reproduce the LS from only the host vari- 
ables. We then build up the Jacobian by taking only 
host variables using both wraparound methods. The 
results are shown in figure ^.b where again the cir- 
cles correspond to square wraparound and the crosses 
to horizontal wraparound. We only plot the first half 
of the spectrum, the second half of the spectrum differs 
considerably for both methods (host-parasites variables 
and only host variables) since the small Lyapunov expo- 
nents are more sensitive to the loss of information con- 
tained in the parasite variables. On the other hand, the 
first half of the spectrum is quite similar independently 
of the choice of host-parasite or only host variables. As 
we can see in figure |l3|.b, horizontal wraparound seems 
to converge to a different curve than square wraparound 
for sub-system sizes N s < 20 (see aligned crosses in the 
lower part of the spectrum). Since we are only taking 
the host population (20 x 20), when N s > 20 the hori- 
zontal wraparound has finished a complete wrap and it 
starts to pick up the neighbours in adjacent rows and 
thus the rescaled spectrum begins to converge closer to 
the square wraparound. 



FIG. 13. First half of the rescaled Lyapunov spectrum 
for a host-parasitoid system in a two-dimensional lattice for 
sub-system sizes N 3 = 1, . . . , 40. a) Using both the host and 
parasite variables and b) using only the hosts when building 
up the Jacobian. The circles (crosses) correspond to the 
square (horizontal) wraparound. 



While the choice of wraparound method is not cru- 
cial for interleaving, figure ^3[a shows that it leads 
to significant differences in rescaling behaviour. Note 
that the LS for the whole system N s = N = 800 is 
not depicted since it would take an enormous amount 
of time to compute. In figure [l^.a we depict the 
rescaled LS for sub-system sizes N s = 1, . . . , 40 for both 
wraparound methods (square wraparound with circles 
and horizontal wraparound with crosses). As for the 
two-dimensional lattice of coupled logistic maps, hori- 
zontal wraparound converges to a different curve than 



IV. CONCLUSIONS 



When studying high dimensional extended dynamical 
systems in a spatio-temporal chaotic regime it is pos- 
sible to rescale the sub-system Lyapunov spectrum to 
obtain the original Lyapunov spectrum. In this ther- 
modynamic limit, a sub-system of comparatively small 
size N s contains a sufficient amount of information to re- 
construct the Lyapunov spectrum of the whole system. 
Usually, when coupling different sub-systems in a lat- 
tice one chooses a coupling with a finite neighbourhood 
(localized coupling) or at least with decreasing effect 
for further away neighbours. In the context of discrete 
spatio-temporal systems, this restriction on the choice 
of coupling causes the Jacobian of the dynamics to be a 
banded (or quasi-banded) matrix. In the limit of only 
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nearest neighbours interaction in a onc-dimcnsional lat- 
tice, the Jacobian is a tridiagonal matrix. If one con- 
siders the homogeneous evolution under this dynamical 
system, the Lyapunov spectrum of sub-Jacobian matri- 
ces will inherit the resc aling and interleaving proper- 



ties described in section II A. The evidence presented 



in this paper shows that the new rescaling method of 
the sub-system Lyapunov spectrum gives a much bet- 
ter fit than the conventional rescaling N/N s for one- 
dimensional lattices. 

We have also observed interleaving of the Lyapunov 
spectra for consecutive sub-system sizes. We showed 
that for two-dimensional lattices the rescaling and inter- 
leaving are still valid. However, the choice of variables 
used to build up the sub-Jacobian matrices appears to 
be crucial to achieve good rescaling properties. In par- 
ticular one has to choose an ordering method of the 
system variables that mimics the propagation of infor- 
mation in the particular lattice topology of the system. 
In two dimensions we showed that choosing the system 
variables in 'concentric' sub-squares gave a much better 
rescaled Lyapunov spectrum than by choosing them in 
a row or column-wise fashion. Generalizing this idea to 
higher-dimensional lattices one should take the system 
variables by filling up 'concentric' hyper-cubes. 

Another point to take into account when choosing 
the system variables in high-dimensional lattices is the 
anisotropy of the coupling. The two-dimensional sys- 
tems studied here have an equal relative contribution 
from all the neighbouring directions (isotropic cou- 
pling). It is possible to choose the coupling in order 
to give more weight to one of the directions (vertical or 
horizontal) and thus the propagation of information to 
be faster in that direction. Therefore, instead of build- 
ing the system variables by 'concentric' squares it should 
be more natural to take rectangles, the ratio of the rect- 
angle sides being related to the ratio of velocity propa- 
gation of disturbances in both directions. 

For a continuous spatio-temporal system a similar re- 
construction may be used by sampling in a grid of a 
sub-system at regular time intervals and by reconstruct- 
ing the Jacobian from time series in the usual manner. 
The same procedure can be applied for a discrete spatio- 
temporal system where the dynamics is not explicitly 
given and the only available dynamic information comes 
from time series taken at several spatial locations. We 
expect that rescaling and interleaving should still be 
observed in these cases. This aspect is currently under 
investigation and will be reported elsewhere 
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